Rcpp::sourceCpp("/home/fra/Documents/GitHub/CommonAtomModel/Functions/CAM/CAM.cpp")
source("/home/fra/Documents/GitHub/CommonAtomModel/Functions/CAM/newCAM.R")
##
## Attaching package: 'reshape'
## The following objects are masked from 'package:plyr':
##
## rename, round_any
## Registered S3 method overwritten by 'pryr':
## method from
## print.bytes Rcpp
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange() masks plyr::arrange()
## x dplyr::combine() masks gridExtra::combine()
## x purrr::compact() masks plyr::compact()
## x purrr::compose() masks pryr::compose()
## x dplyr::count() masks plyr::count()
## x tidyr::expand() masks reshape::expand()
## x dplyr::failwith() masks plyr::failwith()
## x dplyr::filter() masks stats::filter()
## x dplyr::id() masks plyr::id()
## x dplyr::lag() masks stats::lag()
## x dplyr::mutate() masks plyr::mutate()
## x purrr::partial() masks pryr::partial()
## x dplyr::rename() masks reshape::rename(), plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
## Loading required package: coda
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2021 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
6 distribution osservate, 4 gruppi
x0 <- replicate(1, c(rnorm(50,3),rnorm( 150,10,.5)))
x1 <- replicate(1, c(rnorm(100,3),rnorm( 100,10,.5)))
x2 <- replicate(3, c(rnorm(100,-5),rnorm(100,10)))
x3 <- c(rnorm(100,3),rnorm( 100,10,.5),rnorm( 100,-5,.5),rnorm( 100,0,.5))
yd <- c(x0,x1,x2,x3)
yg <- rep(1:6,c(rep(200,5),400))
plot(yd,col=yg)
m1 <- CAM(y_obser = yd,
y_group = yg,
K0 = 20,
L0 = 20,
prior = list(
# hyperparameters NIG
m0= mean(yd),
k0=1/var(yd),
a0=3, b0=1,
# hyperparameters alpha and beta
a_alpha=3, b_alpha =3,
a_beta =3, b_beta = 3),
nsim = 20000,
burn_in = 20000,
thinning = 1,
verbose = F,
fixedAB = F,
kappa=0.5,
cheap=F,
seed=123456,sigma.prop.beta = 1
)
plot(ts(m1$A_DP))
acf(m1$A_DP)
pacf(m1$A_DP)
plot(ts(m1$B_DP))
acf(m1$B_DP)
pacf(m1$B_DP)
r1=mcclust::comp.psm(as.matrix(t(map_dfc(m1$Z_j, ~.x) )))
## New names:
## * NA -> ...1
## * NA -> ...2
## * NA -> ...3
## * NA -> ...4
## * NA -> ...5
## * ...
image(r1)
a1=mcclust.ext::minVI(r1,method = "greedy")
plot(yd,col=rep(a1$cl,table(yg)))
r=mcclust::comp.psm(as.matrix(t(map_dfc(m1$Csi_ij, ~.x) )))
## New names:
## * NA -> ...1
## * NA -> ...2
## * NA -> ...3
## * NA -> ...4
## * NA -> ...5
## * ...
image(r)
a=mcclust.ext::minVI(r)
plot(yd,col=a$cl)
xx <- seq(-10,15,by = .05)
LLL <- list()
LLL <- posterior.densities(m1,howfarback = 1000,T,xx = xx)
## 123456
par(mfrow=c(1,2))
for(hh in 1:6){
hist(yd[yg==hh],breaks = 30,freq=F, main="hist+mcmc+mean")
for(i in 1:1000)
lines(LLL[[hh]][i,]~xx,col=1)
lines(colMeans(LLL[[hh]])~xx,col=2,lwd=3)
hist(yd[yg==hh],breaks = 30,freq=F, main="hist+mean")
lines(colMeans(LLL[[hh]])~xx,col=2,lwd=3)
}
# sim <- 500
# ind <- 500
# ygind <- yg[[ind]]
#
# plot(yd)
# points(yd[ind]~c(ind),col=2,pch=21,bg=2)
#
# mat <- matrix(NA,nsim,2)
# for(sim in 2:nsim){
#
# mat[sim,1] <- m1$OMEGA[[sim]][ m1$Csi_ij[[sim]][ind]
# ,
# m1$Z_j[[sim]] [ygind]]
# mat[sim,2] <- m1$THETA[[sim]][ m1$Csi_ij[[sim]][ind]
# ,
# 1]
# }
#
# plot(ts(mat))